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Abstract 


We investigate several distribution free dependenee deteetion proeedures, mainly based 
on bootstrap prineiples and their approximation properties. Thanks to this study, we 
introduee a new distribution free Unitary Events (UE) method, named Permutation UE, 
whieh consists in a multiple testing procedure based on permutation and delayed co¬ 
incidence count. Each involved single test of this procedure achieves the prescribed 
level, so that the corresponding multiple testing procedure controls the Ealse Discovery 
Rate (EDR), and this with as few assumptions as possible on the underneath distribu¬ 
tion. Some simulations show that this method outperforms the trial-shuffling and the 
MTGAUE method in terms of single levels and EDR, for a comparable amount of false 
negatives. Application on real data is also provided. 


1 Introduction 


The eventual time dependence either between cerebral areas or between neurons, and 
in particular the synchrony phenomenon, has been vastly debated and investigated as a 


potential element of the neuronal code (Singer 19931. To detect such a phenomenon 


at the microscopic level, multielectrodes are usually used to record the nearby electri¬ 
cal activity. After pretreatment, the time occurrences of action potentials (spikes) for 
several neurons are therefore available. One of the first steps of analysis is then to un¬ 
derstand whether and how two simultaneously recorded spike trains, corresponding to 
two different neurons, are dependent or not. 
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Several methods have been used to deteet synehrony (Perkel et al. 1967; Aertsen 
et al411989|). Among the most popular ones, the UE method, due to Griin and eol- 


laborators (Griin 1996 Griin et al. 20101, has been applied in the last deeade on a 


vast amount of real data (see for instance (Kilavik et al. 2009) and references therein). 


Two of its main features are at the root of its popularity: UE is not only able to give a 
precise location in time of the dependence periods, but also to quantify the degree of 
dependence by providing p-values for the independence tests. 

Erom the initial method, substantial upgrades have been developed: 

(i) In the original UE method, the point processes modelling the data are binned and 
clipped at a rough level, so that the final data have a quite low dimension (around a few 


hundreds per spike train). However, it is proved in (Griin et al. 1999) that the binned 
coincidence count as a result of this preprocessing may induce a loss in synchrony 


detection of about 60% in certain cases. The idea of (Griin et al. 1999) was therefore 


to keep the data at the initial resolution level despite its high dimension, but to define 
the notion of multiple shift (MS) coincidence count, nicely condensing the dependence 
feature that neurobiologists want to analyze without any loss in synchrony detection. 

(ii) The original UE method assesses p-values by assuming that the underlying point 
processes are Poisson (or Bernoulli) processes. As there is still no commonly validated 


and accepted model for spike trains (see for instance Pouzat & Chaffiol (2009) for 


a thorough data analysis), several surrogate data methods have been proposed (Eouis 


et al.[ |2010[ ). In particular, trial-shuffling methods (Pipa & Grun[ |2003t |Pipa et ah 


2003) allow to assess p-values based on the fact that i.i.d. trials are available, using 
bootstrap paradigm and without making any assumption on the underlying point pro- 


3 









































cesses distribution. However, up to our knowledge, surrogate data methods are always 
based on binned coincidence count (see section for a precise definition) whose low 


complexity combined to parallel programming (Denker et al. 20101 make algorithms 


usable in practice. 

(iii) The original UE method is based on two main statistical approximations. First, 
it involves the underlying intensities (or firing rates) of the Poisson processes, which 
are unknown in practice, and so replaced by estimates. However, this plug-in proce¬ 
dure is not taken into account in the p-values computations. Then, the detection of 
dependence through time is done by applying several tests at once without correcting 


for the multiplicity of the tests. In the recent work of (Tuleau-Malot et al. 2014), a 


multiple testing procedure based on a Gaussian approximation of the Unitary Events 
(MTGAUE) corrects those two facts, moreover including a generalization of the notion 
of MS coincidence count, namely the delayed coincidence count, which does not suffer 
from any loss in synchrony detection. But MTGAUE is still based on the assumption 
that the underlying point processes are Poisson. 

Our aim is here to go further by proposing a new delayed coincidence count-based 
multiple testing procedure, which does not need any binning preprocessing of the data 


as in (Tuleau-Malot et al. 20141, but which does not assume any model on the underly¬ 


ing point processes anymore. This procedure combines a permutation approach in the 
line of ( IHoeffdiri^ [1952 Romano 1989 Romano & Wolf 2005), with the multiple 
testing procedure of ( [Benjamini & Hochberg[|1995| ). 

To do so, we first propose a fast algorithm to compute the delayed coincidence 
count, with a computational cost equivalent to the one of the binned coincidence count. 
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Next we study several distribution free tests, most of them based on bootstrap ap¬ 
proaches, as the trial-shuffling, or the permutation approach. We finally propose a 
complete multiple testing algorithm which satisfies similar properties as existing UE 
methods, but without sharing any of the previous drawbacks. Some simulations and 
applications to real data complete this study. 

In all the sequel, and denote two point processes modelling the spike trains 
of two simultaneously recorded neurons and X represents the couple (X^,X^). The 
abbreviation ’’i.i.d.” stands for independent and identically distributed. In this sense, by 
assuming that n independent trials are observed, the observation is modeled by an i.i.d. 
sample of size n of couples from the same distribution as X, meaning n i.i.d. copies 
Xi, ...,Xn of X. This sample is denoted in the sequel by = (Xi, ...,X„). The 
corresponding probability and expectation are respectively denoted by P and E. 

Since the independence between X^ and X^ is the main focus of the present work, 
the following notation is useful: X^ denotes a couple (X^’^, X^’^) such that X^’^ 
(resp. X^’^) has the same distribution as X^ (resp. X^), but X^’^ is independent of 
X^’^. In particular, the couple X^ has the same marginals as the couple X. Moreover, 

denotes an i.i.d. sample of size n from the distribution of X^, and Px and Ex are 
the corresponding probability and expectation. 

Note in particular that if the two observed neurons indeed behave independently, 
then the observed sample X„ has the same distribution as X^. 

The notation lx stands for a function whose value is 1 if the event A holds and 0 
otherwise. 

Finally, for any point process X^ (j = 1, 2), dNxi stands for its associated point 
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measure, defined by: 


f{u)dNxj {u) = ^ /(T) for all measurable funetion /, 


T&Xo 


and for any interval I, Nxo (/) denotes the number of points of observed in /. 


2 Binned and delayed coincidence counts 

Beeause of the way neurons transmit information through action potentials, it is com¬ 
monly admitted that the dependence between the spike trains of two neurons is due to 
temporal correlations between spikes produced by both neurons. Informally, a coinci¬ 
dence occurs when two spikes (one from each neuron) appear with a delay less than 
a fixed 5 (of the order of a few milliseconds). Several coincidence count functions 
have been defined in the neuroscience literature, and among them the classical binned 


coincidence count, used for instance in ([Pipa & Grun[ |200^ |Pipa et al.[|2003[). 


Definition 1 The binned coincidence count between point processes and X"^ on the 
interval [a, b] with b — a = M5 for an integer M > 2 and a fixed delay 6 > 0 is given 
by 

M 

£=l 

where Ii is the ith bin of length 5, i.e. [a -f (f — 1)5, a + £6) and 




1 if there is at least one point of XI in the ith bin, 
0 if there is no point of X^ in the ith bin. 


More informally, the binned coincidence count is the number of bins that contain at least 
one spike of each spike trains. The binned coincidence count computation algorithm is 
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usually performed on already binned data. Therefore, given two sequenees of 0 and 1 of 


length {b — a)6 ^,2{b — a)6 ^ operations are needed to eompute the binned eoineidenee 
count, without counting the number of operations needed for the binning preprocessing. 


The more recent notion of delayed coincidence count, introduced in (Tuleau-Malot 


et al. 20141, is a generalization of the multiple-shift coincidence count, defined in (Griin 


et al. 19991 for discretized point processes, to non necessarily discretized point pro¬ 


cesses. 


Definition 2 The delayed coincidence count between point processes and on 
the interval [a, 6] is given by 

^cornc^X\X^)= [ [ l|„_,|<,diVxi(M)diVx2(n), 

J a J a 

More informally, X^) is the number of couples of spikes (one spike from Xi 

and one from X 2 ) appearing in [a, b] with delay at most equal to 6. 

The delayed coincidence count c := X^) can be computed using the 

following algorithm. 
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Delayed coincidence count algorithm 


Given two sequences xi and X 2 of ordered points with respective lengths ni = 
NxT-{[a,b]) and n 2 = Nx^da^b]), representing the observations of two point pro¬ 


cesses and 

- Initialize j = I and c = 0. 

- For i = I,Til, 

1. Assign xiow = xi[i] - 6. 

2. While j < n 2 and X 2 [j] < xiow, 3=3 + 1- 

3. If J > n 2 , stop. 

4. Else (here necessarily, X 2 [j] > xiow), 

4.a Assign Xup = Xi [i] + 5 and k = j. 

4.b While k <n 2 and X 2 [k] < Xup, c = c + 1 and k = k + 1. 


It is easy to see that the complexity of this algorithm is not governed only by the 
lengths Til and n 2 but also by the random numbers of points in intervals of length 25 
for step 4.b. More precisely, it is bounded by 3ni (for steps 1, 3 and 4.a), n 2 (for 
all steps 2 on all points of xi) and 2ni times the number of points of X 2 in a segment 
(namely [xiow, Xup\) of length 25 (for step 4.b). In average, if and are for instance 
independent homogeneous Poisson processes of respective intensities Ai and A 2 , at most 
3Ai(6 — a) + A 2 (& — a) + 2Ai(6 — a) ( 25 X 2 ) operations are required. So, for typical 
parameters ((& — a) = 0.1s, 5 = 0.005s, Ai = A 2 = 50Hz), 40 operations in average 
are required to compute the binned coincidence count, against 25 operations for the 
delayed coincidence count. Both algorithms are therefore linear in (b — a) with a slight 
advantage for the delayed coincidence count which exploits the sparsity of the spike 
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trains, in the usual range of parameters. Notiee that all surrogate data methods (see 


(Louis et al. 20101) eould in prineiple be applied on this new notion of ooineidenee, at 


least when only two neurons are involved. 


3 Some distribution free independence tests 


Given the observation of a n sample = (Xi,..., X„) eorresponding to n different 
trials, the aim is to test: 


and are independent on [a, 6]’ 


against 


{Hi) ”X^ and X^ are not independent on [a, 6]”. 


All existing UE methods are based on a statistie equal to the total number of eoinei- 


dences: 


C = C(X„) = 5^V>(A-‘,At), 


2 = 1 


where tp generically denotes either or or other eoincidence eount func¬ 


tions that practitioners would like to use (see (Albert et al. 2014) for other choices). 


To underline what is observable or not, when C is computed on the observation of 
X„, it is denoted by the total number of observed coincidences. 

In the following, several of these UE methods or testing procedures are described, 
which all rely on the same paradigm: ’’reject (Hq) when is significantly different 
from what is expected under {Hq)”. More precisely, the independence {Hq) is rejected 
and the dependence is detected when a quantity, based on the difference between the 
observed coincidence count and what is expected under (Hq), is smaller or larger than 
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some critical values. Those critical values are obtained in various ways, each of them 
being peculiar to each method. Note that the following procedures could be applied to 
any chosen coincidence count function 99 , though the implemented procedures of the 
present simulations and data analysis only use the delayed coincidence count, that is 

(p = 


3.1 A naive approach and the centering issue 

As noticed above, the only question, when considering UE methods, is how to construct 
the critical values. 

If the values of the expectation and the variance of C under (Hq), that is 
Co = Ex(C) and vq = Ex ((C - cq)^) , 


are precisely known, then the classical Central Limit Theorem gives under indepen¬ 
dence that 


C(X;^)-co ^ 

n^o 


Ar(0,l). 


Then, given a in (0,1), the test which consists in rejecting {Hq) when — cq)/a/uq 
is larger than Zi_a, the 1 — a quantile of a standard Gaussian distribution, is asymptot¬ 
ically of level a. It means that, for this test, the probability of rejecting independence, 
whereas independence holds, is asymptotically (in n, the number of trials) smaller than 
the prescribed a. 

In the present point processes framework, strong distribution assumptions, for which 
the values of Cq and Vq are precisely known, are unrealistic. Even if the spike trains are 


assumed to be homogeneous Poisson processes as in (Tuleau-Malot et aL| 20141, those 
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quantities depend, through some formulas, on the unknown firing rates that have to be 


estimated and plugged into these preeise formulas. It has been shown that this modi¬ 
fies the asymptotie varianee shape, and tests under Poisson assumptions with unknown 


firing rates have been developed in (Tuleau-Malot et al. 2014). 


Since Poisson assumptions are quite doubtful on real data (Pouzat & Chaffiol, 20091, 


the aim of the present work is to go further by not assuming any Poisson or other model 
assumptions for the spike trains. In this sense, the aim is to develop ’’distribution free” 
methods that are completely agnostic with respect to the underlying distribution of the 
spike trains. In this case, a preliminary step is to estimate cq, only using the sample 
without any distribution assumption. Note that 


Co = Ejl 


n 

E 

2=1 






and that for i ^ i' ,d& Xi is always assumed to be independent of Xj/, 


E [tp(XlXl)\ = Ei , 


( 1 ) 


Therefore, Cq can always be estimated in a distribution free manner by 


C„(x„) = —y]4,(A‘,A.?) 


i^i' 


Hence a reasonable test statistic would be based on the difference: 


U = U(X„) = C(XO-Co(X„), 

its observed version being denoted by Here, U(X„) is not an empirical mean, 

but a 17-statistic, so it does not satisfy the classical Central Limit Theorem. Hence, 
its limit distribution under (Hq) is not as straightforward as usual. Nevertheless, some 
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asymptotic theorems, proved in (Albert et al. 2014), show that under mild eonditions 
(always satisfied in praetiee in the present eases) the following eonvergenee result holds: 

U(X;[) 


Z(X;^) = 


V^<t(X;^) 




( 2 ) 


where 




n{n — l){n — 2) 


h{X„X,)h{Xi,Xk), 


i,j,k all different 


with 


y) = n ) + ^iy )-^iy ) 


2 L 

As above, denoting by the quantity Z eomputed on the observed sample, (|^ 
implies that for some fixed a in (0,1), the test that eonsists in rejeeting (Hq) when 
2 ^obs y zi-a, is asymptotieally of level a. 

The approximation properties of ([^ are illustrated on Figure 

Clearly, one ean see that the distribution approximation is good when n is large 
(n = 200) as expected, but not so eonvincing for small values of n (n = 20, or even 
n = 50), partieularly in the tail parts of the distributions. However, as it is espeeially 
the tails of the distributions that are involved in the test through the quantile Zi^a, one 
can wonder, by looking at Figure if it may perform reasonably well in practice with 
a usual number of a few tens of trials. 

Furthermore, looking informally at Equation ([^, readers eould think of two approx¬ 
imations that eould be roughly formulated in the following way: 


U(X: 


J-N 


Ar(0,n(T^(X;^)) , 


(3) 


and 


C(X;^) S V(C„(Xi-),na^(X^) 


-L^ 




(4) 
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n = 20 


n = 50 


n = 200 





Figure 1: Gaussian approximation of the distribution of Z. In plain blaek, eumulative 
distribution funetion (e.d.f.) of Z under (Hq), that is of Z^ = Z(X;J-) obtained with 2000 
simulations of X^, for n = 20, 50 or 200 trials of two independent Poisson proeesses 
of firing rate 30Hz, on a window of length 0.1s with 6 = 0.01s. The dashed line 
eorresponds to the standard Gaussian e.d.f. 

This is illustrated on Figure]^ 

Looking at the first line of Figure one ean see that the approximation formulated 
in Q is aetually eoneeivable for large values of n. Note that in praetiee, one eannot 
have aeeess to (5-^(X^) and it has to be replaeed by (j^(X„), meaning that it is eomputed 
with the observed sample. This does not ehange anything under (Hq) sinee X„ is in this 
ease distributed as X;^. But this is a partieularly important stieking point if (Hq) is not 
satisfied as one ean see on the third line of Figure the distribution of U(X;J-) does 
not look like a eentered Gaussian distribution of varianee n(T^(X„), when X„ does not 
satisfy (Ho). 
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Figure 2: Other Gaussian distribution approximations. Two first lines: c.d.f. of U and 
C under (Hq), obtained as in FigureThese c.d.f. are respectively compared with the 
Gaussian c.d.f. with mean 0 and standard deviation y/na{Xn), and the Gaussian c.d.f. 
with mean Co(X„) and standard deviation ^/n^^{Xn), for five different simulations of 
X„ under (Ho). Third line: c.d.f. of U under (Hq) computed as above, compared 
with the centered Gaussian c.d.f. with standard deviation ^/na{Xn), for five different 
simulations of X„ under (Hi) (same marginals as in the first two lines but = X^). 
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More importantly, the seeond line of Figure]^ shows that the approximation formu¬ 
lated in Q is in faet misleading. To understand why, one needs to take into aeeount the 
two following points. 

(i) Co(X;J-) moves around its expeetation cq (whieh is also the expectation of C(X;^)) 
with realizations of X^. These fluctuations have an order of magnitude of ^/n and are 
therefore perfectly observable on the distribution of C(X;J-) whose variance is also of 
order ^/n. 

(ii) n(T^(X;J-) estimates the variance of U(X^) and not the one of C(X;J-) or Co(X^). 
This explains why not only the mean but also the variance are badly estimated in the 
second line of Figure]^ Both randomness (the one of C(X^) and the one of Co(X;J-)) 
have to be taken into account to estimate the variance of U (X^). 

As a conclusion of this first naive approach, the test of purely asymptotic nature, 
which consists in rejecting (Hq) when > zi^a may work for n large enough, as 
the variance is here computed by considering the correctly recentered statistic U, and 
this even if the behavior of the statistic under (Hi) is not clear. But an ad hoc and more 
naive test statistic, based on an estimation of the variance of C directly and without 
taking into account the fact that the centering term Co(X„) is also random, would not 
lead to a meaningful test. 


3.2 The bootstrap approaches 


It is well known (Gine 19971 that tests of purely asymptotic nature as the one presented 


above are less accurate for small values of n than more involved procedures. In this arti¬ 


cle, the focus is on bootstrap/resampling procedures that are usually known to improve 
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the performance from moderate to large sample sizes. Three main procedures are in¬ 


vestigated: the trial-shuffling introduced in (Pipa & Griin 2003; Pipa et al. 20031, the 


classical full bootstrap of independence and the permutation approach (Romano 19891. 


The main common paradigm of these three methods, as described in the sequel, 
is that starting from an observation of the sample X„, they randomly generate via a 
computer another sample X„, whose distribution should be close to the distribution of 
X;J- (see also Figure]^. 

^ Trial-shuffling ^ 


X„ = Xf = ((A'.W,1,, .... 

where the (k), {k)ys are n i.i.d. couples drawn uniformly at random in 

{{hj) / * = ^ j}. 

In particular, the corresponding bootstrapped coincidence count is 


CT5 ^ C(X™) := J . 

k=l 


This algorithm seems natural with respect to ([T]) because it avoids the diagonal terms 
of the square {(?, j) / i = 1, j = 1, Hence as a result. 


E(C^^) = Co = Ex(C). 
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r Classical full bootstrap ^ 

X,. = x; = 

where the n couples {i* {k), j* {k)) are i.i.d. and where i*{k) and j*{k) are drawn 
uniformly and independently at random in {1,n}. 

In particular, the corresponding bootstrapped coincidence count is 


n 

C- = C(X;) 


k=l 


Note that this algorithm draws uniformly at random in the square {{i,j) / i = 
j = and therefore does not avoid the diagonal terms. The idea behind 

this algorithm is to mimic the independence under (Hq) of Xl and by drawing the 
indexes i*{k) and j*{k) independently. However 


E(C*) =n 


n n 


Hence under (Hq), Ex(C*) = Cq but, under (Hi), E(C*) and cq are only asymptotically 
equivalent. 

Permutation ^ 


where is a permutation drawn uniformly at random in the group of permutations 
&n of the set of indexes n}. 

In particular, the corresponding bootstrapped coincidence count is 


C* = C(x;~) 


2=1 
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The idea is to use permutations to avoid pieking twiee the same spike train of the 


same trial. In partieular under (Hq), the sum in C* is still a sum of independent vari¬ 
ables, whieh is not the ease in both of the previous algorithms. However, under (Hi), 
the behavior is not as limpid. As for the full bootstrap. 


E(C*) = n 


1 in _ 1 

+ - Ej_{^{X\X^)) 

n n 


Henee under {Hi), E(C*) and cq are only asymptotieally equivalent. 


To eompare those three bootstrap/resampling algorithms, the first thing to wonder 
is whether, at least under {Hq), the introdueed extra randomness has not impaeted the 
distribution. More preeisely, as stated above, all three proeedures satisfy 

Ex(C(X„)) = Ex(C(X^)) = Co, 

but is the full distribution of C(X„) the same as the one of C(X;J-)? 

The first line of Figure shows as expected that the permutation does not change 
the distribution of X^, since, as said above, no spike train is picked twice. However, 
clearly the trial-shuffling and the full bootstrap have not the same property, even if the 
distributions are quite close. 

Nevertheless, this is not completely convincing. Indeed, the main advantage of 
bootstrap procedures is to be able for one current observation of X„ to generate several 
realizations of X„ to obtain not the unconditional distribution of C(X„) but the condi¬ 
tional distribution of C(X„) given X„. Figure gives a more visual representation of 
the difference between conditional and unconditional distributions. 
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Conditional Unconditional 


Trial-Shuffling 


Full Bootstrap 


Permutation 



Figure 3: The uneonditional distribution and eonditional distributions of C under {Hq). 
C.d.f. of C(X;^) and (for the first line) of = C(X^^), of C* = C(X;) and of 
C* = C(XJJ”) obtained from 10000 simulations of n = 20 trials of two independent 
Poisson proeesses of firing rate 30Hz on a window of length 0.1s with 5 = 0.01s. On the 
seeond line, in addition to the e.d.f. of C (X^), five observations of X„ = X;J- have been 
simulated in the same set-up and given these observations, the eonditional e.d.f. have 
been approximated by simulating 10000 times the extra-randomness eorresponding to 
X„. For the trial-shuffling, in addition to this approximate Monte-Carlo method (MC), 


the exaet eonditional e.d.f. has been obtained thanks to the algorithm of (Pipa et al 


2003). 
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Figure 4: Schematic view of the bootstrap paradigm and the difference between un¬ 
conditional distribution and conditional distribution given the observation. 

In particular, in a bootstrap testing procedure, the critical values are quantiles of 
the conditional distribution of the bootstrapped test statistic given the observation of 
X„ and not the quantiles of the unconditional distribution. Hence, to see whether boot¬ 
strapped critical values associated to C(X„) are reasonable, the conditional distribution 
of C(X„) given X„ has to be compared with the distribution of C(X^), and this whether 
X„ satisfies (Hq) or not. 

^In fact, the quantiles are usually approximated by a Monte-Carlo method, since one has access to 
only a huge but still finite number of realizations of X„ given X„ in practice. 
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The second line of Figure shows what happens for the approximated conditional 
distribution of C(X„) given X„ under (Hq) in the three considered bootstrap approaches. 
Surprisingly none of these three conditional distributions seems to fit the distribution of 
C(X^). One may eventually think that this is due to the Monte-Carlo approximation 
of the conditional distributions, but for the trial-shuffling approach, Pipa and Griin de¬ 


veloped an algorithm for exact computation of the conditional distribution (Pipa et al 


2003): both Monte-Carlo and exact conditional distribution are so close that it is diffi¬ 


cult to make any difference between them. Hence there should be another explanation. 
In fact, the curves on the second line of Figure are similar to the ones on the second 
line of Figure]^ In both set-ups, one wonders if the distribution of C(X^) can or can¬ 
not be approximated by a distribution depending on the observation of X„: a very basic 
Gaussian distribution for Figure]^ and a more intricate distribution using the bootstrap 
paradigm for Figure Nevertheless both distributions are too widely spread around 
the aim which is the distribution of C(X;J-). Since the explanation for Figure was a 
centering defect that can be corrected by considering U, one of the possible explanation 
here is a centering defect for the bootstrap procedures too. 


3.3 Which centering for which bootstrap ? 


All the bootstrap approaches that have been proved to work from a mathematical point 


of view are based on centered quantities (Gine, 19971. In particular, the precursor work 


of Bickel and Freedman (Bickel & Freedman 19811 on the bootstrap of the mean can 


be heuristically explained as follows. 

Given a n sample of i.i.d. real random variables = (Ti,..., Yn) with mean m 
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and a corresponding bootstrap sample Y* , it is not possible to estimate the distribution 
of the empirieal mean V = (1/n) direetly. However one ean estimate the 

centered distribution, i.e. the distribution ofY — m = Y — E(y). To do so, it is 
suffieient to replaee ’’empirical mean” by ’’empirieal bootstrap mean” and ’’expeetation” 
by ’’conditional expectation”. More explieitly, denoting by Y* the empirical mean of the 
bootstrap sample Y* , the distribution of F — E(F) is approximated by the eonditional 
distribution given Y„ of Y* — E(F* |Y„). 

Transposed in our framework, this paradigm would mean that one ean obtain a 
good fit of the distribution of (l/n)(C(X;J-) — cq) by the eonditional distribution of 
(l/n)(C(X„) — E(C(X„)|X,^)) given X„. But as seen above, construeting a test with 
test statistie (1/?7 ,)(C(X„) — cq) is impossible in a full distribution free eontext where 
the value of Cq is unknown. 

Therefore one needs to find a quantity that is eompletely observable but whose mean 


is still null under {Hq). The statistie U introdueed in Section 3.1 is suitable from this 
point of view. What one needs to check is whether the distribution of U(X„) under 
{Ho), that is of U(X^) (whieh has zero mean), is well approximated by the distribution 

ofu(x„) -e(u (x„) |X,). 

For the trial-shuffling, sinee 


u(xr) = y;v>(.vw,.A'| 


k=l 


’(fc) 


n — 1 




k^k' 


TS[k') I ) 


one ean easily see that beeause the eouple [i'^^{k),j'^^{k)) is drawn uniformly at ran- 
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dom in the set of the (i, j)’s sueh that i ^ j (set of eardinality n{n — 1)), 

E(u(xr)|x„) = 

Co (X„) - C (X„) 
n 

U(X„) 

n 

Henee the eorreet bootstrap statistic is 

_ jj = u (X™) + ^ ■ 

However similar computations show that the full bootstrap and the permutation satisfy 


E(U(X:)|XO=E(U (X^)|X„) =0, 

so U(X;)andU(Xj") can be used directly. 

Figure shows the quality of approximation of the distribution of U (X;J-) by the 
conditional distribution given the observation of either U* = U (X*)orU* = U (X^"). 
The approximation is accurate under (Hq) but it is actually also accurate even if the ob¬ 
served sample is simulated under (Hi), which is in complete accordance with the math¬ 


ematical results of consistence in Wasserstein distance proved in (Albert et al. 2014). 
The approximation is just as accurate for the recentered statistic -f 

Note that the difference between the conditional c.d.f. of and the one of is 
particularly visible under (Hi) when = X^. Hence, as explained by the compu¬ 
tations above, in a trial-shuffling approach, the recentered version leads to the correct 
bootstrap distribution. Note finally that this corroborates the previous intuition: the rea¬ 
son why the approximation works for U and not for C is exactly the same as for the 
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Trial-Shuffling 


Full Bootstrap 


Permutation 



Figure 5: Conditional distribution of U (or its recentered version for the 
Trial-Shuffling) given X„. Cumulative distribution functions ofU^ = U (X;^) in black, 
obtained by simulation as in Figure For the first line, under (Hq), five observations 
of X;^ in the same set-up have been fixed and given these observations, the conditional 
c.d.f. of = U (X™), of + U'^^Vn, of U* = U (X;) and of U* = 

U (X^”) have been obtained as in Figure]^ For the second line, five observations of 
X„, simulated under (Hi) with marginals equal to the ones of the first line but satisfying 
have been simulated and conditional c.d.f. are obtained in the same way as 

above. 

naive approach of Figure The centering is indeed random (here it can be viewed as 
E(C(X„) |X„)) and one needs to take it into account to have a correct approximation. 


24 



























































Finally an extra simplification holds in the permutation case, which may seem very 
surprising. One can easily rewrite on the one hand, 

u (X») = (i - c (X„) - E ^ E) 

and, on the other hand, for the permutation sample 

u (x"“) = (l - c «") - E*’ W-E) ■ 

^ ^ ij 

Note that the sum ^ (p (X/, X|) is invariant by the action of the permutation. Hence 
if M* denotes the quantile of order t of the conditional distribution of U given 

X„ and if Cj denotes the quantile of order t of the conditional distribution of C (Xjj'") 
given X„, this very simple relationship holds 

^ ^ hj 

Hence the test that rejects (Hq) when U (X„) > is exactly the one that rejects 
(Hq) when C (X„) > c\_^. Therefore despite the fact that the conditional distribution 
of C (Xjj") is not close at all to the one of C (X;^), the test based on C works, because 
it is equivalent to the test based on U, for which the approximation of the conditional 
distribution works. Note however that this phenomenon happens only in the permuta¬ 
tion approach, but not in the trial-shuffling or the full bootstrap approaches. 

3.4 Practical testing procedures and j9-values 

From the considerations given above, five different tests may be investigated, the first 
one based on a purely asymptotic approach, and the four other ones based on bootstrap 
approaches, with critical values approximated through a Monte-Carlo method. For each 
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test, the eorresponding p-values (i.e. the values of a for whieh the test passes from 
aeoeptanee to rejeetion) are given. 

The naive test (N) It eonsists in rejeeting (Hq) when 

The corresponding p-value is given by: 

1 _ $ , 

where $ is the c.d.f. of a standard Gaussian distribution. 


The Trial-Shuffling test, version C (TSC) It consists in rejecting (Hq) when 


Qiobs 


— G-o 




where is the empirical quantile of order (1 — a) of the conditional distribution 
of given X„. This empirical quantile is estimated over B {B = 10000 usually) 
realizations C™,..., given the observed sample X„. The corresponding p-value is 
given by: 


B 


B 


i=l 




Despite the problems underlined in Section 3.3 we kept this test in the present study 


since it corresponds to the one programmed in (Pipa & Griin 20031 and since this 


bootstrap procedure is usually the one applied by neuroscientists. 


The Trial-Shuffling test, version recentered U (TSU) It consists in rejecting (Hq) 
when 
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where is the empirieal quantile of order (1 — a) of the eonditional distribution 
of (the eorreetly reeentered statistie) given X„. This empirieal quantile and the 
corresponding p-value are obtained in a similar way as for the above (TSC). 

The Full Bootstrap test, version U (FBU) It consists in rejecting {Ho) when 

where ul_^ is the empirical quantile of order (1 — a) of the conditional distribution of 
U* given X„. This empirical quantile and the corresponding p-value are obtained in a 
similar way as for the above (TSC). 


The permutation test (P) The reader may think that it should consist in rejecting 
(Hq) when 

f-^obs \ 


where is the empirical quantile of order (1 — a) of the conditional distribution of 
C* given X„. But the test by permutation is in fact directly defined by its p-value, which 
is slightly different here, equal to: 


B 






i=l 


B + 1 

The permutation test then consists in rejecting (Hq) when this p-value is less than a. 
Indeed, such a permutation test, with such a slightly different version of p-value, has 
been proved to be exactly of level a, whatever B ( jRomano & Wolf 20051. 

Note however that such a slight correction does not work for full bootstrap or trial¬ 


shuffling approaches, where the tests are only guaranteed to be asymptotically of level 


a. 
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Figure 6: Distribution of the p-values for the different tests. C.d.f. under both (Hq) 
and (Hi) of the p-values for the five tests: naive (N), Trial-Shuffling version C (TSC), 
Trial-Shuffling version U (TSU), Full Bootstrap version U (FBU), and Permutation (P). 
Under (Hq), the c.d.f. are obtained by simulations done as in Figure]^ the c.d.f. are then 
plotted only for small p-values (< 0.1). Under (Hi), the couple (X^, is constructed 


by injection (Griin et al. 2010 Tuleau-Malotetal. 2014l,i.e. as 


where (X^, X^) are two independent Poisson processes of firing rate 27 Hz on a window 
of length 0.1s and where is a common Poisson process of firing rate 3Hz; once 
again, 20 i.i.d. trials are simulated 10000 times to obtain the corresponding c.d.f. with 


5 = 0.01s. 
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Saying that a test rejects at level a is exactly equivalent to saying that its p-value is 
less than a. If a test is of level a for any a in (0,1), the c.d.f. of its p-values should 
therefore be smaller that the one of a uniform variable (i.e. the diagonal). Between 
several tests with this guarantee, the less conservative one is the one for which the c.d.f 
of its p-values is the closest to the diagonal. The left hand-side of Figure shows the 
c.d.f. under (Hq) of the corresponding p-values for the five considered testing proce¬ 
dures and focuses on small p-values, which are the only ones usually involved in testing, 
to highlight the main differences between the five methods. For the chosen value of n 
(n = 20), the c.d.f. of the (TSU) and (FBU) p-values are almost identical and above 
the diagonal, meaning that the corresponding tests do not guarantee the level. On the 
contrary, the c.d.f. of the (N) and (TSC) p-values are clearly under the diagonal and 
far from it, meaning that the corresponding tests are too conservative. As guaranteed 
by ( [Romano & Wolf [[2005 1 ), the permutation approach guarantees the level of the test: 
the c.d.f. of the (P) p-values is also under the diagonal, but much closer to the diagonal 
than the one of the (N) and (TSC) p-values. 

Furthermore, the behavior of the c.d.f. of the p-values under (Hi) gives an indica¬ 
tion of the power of the test: the highest the c.d.f. of the p-values, the most powerful the 
corresponding test. Hence among the tests that guarantee the level, the permutation test 
(P) is the most powerful one. Note that other simulations in more various cases have 


been performed in (Albert et al. 20141 leading to the same conclusion. 


In the sequel, the focus is therefore on the permutation approach, keeping also the 
trial-shuffling version C approach as a variant of the method developed in ( [Pipa & Grtin 


2003). 
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4 Multiple tests 


4.1 Description of the complete multiple testing algorithm 

To detect precise locations of dependence periods that can be matched to some ex¬ 
perimental or behavioral events, it is classical to consider a family of windows W of 
cardinal K, which is a collection of potentially overlapping intervals [a, b] covering the 


whole interval [0,T] on which trials have been recorded (Griin et al. 1999; Tuleau 


Malot et al. 20141. Then, some independence tests are implemented on each window 


of the collection. Here we propose a complete algorithm which takes into account the 
multiplicity of the tests, and which moreover enables to see if the coincidence count is 


significantly too large or too small on each window as in (Tuleau-Malot et aL| 20141. 
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Permutation UE algorithm 


Fix a real number g in (0,0.5) and an integer B larger than 2. 

- Do in parallel for all window W = [a, 6] in W: 

* Extract the points of the X/’s and Xf's in [a,b]. 

* For all (i, j) in {1,compute aij = (X/, over [a, b] 

by the delayed coincidence count algorithm. 

* Draw at random B i.i.d. permutations Ilj!), 1 < b < i?, and compute a, nb(i). 

* Compute also 

* Return ^ ('l + lcb>co6.) andp^^ = ;^ (l + Ef=i lcb<c°o-)- 


- Perform the procedure of (Benjamini & Hochberg 1995 1 on the set of the above 2K p-values: 

* Sort the p-values < <p(2^). 

* Find k = max{Z / < lq/{2K)}. 

* Return all the (PF, ewY^, for which W is associated with one of the p-values p^^^ for I < k, 

with eiy = 1 if p(^ < so the coincidence count is significantly too large on W, 
and ew = — 1 if Pw — coincidence count is significantly too small on W. 


This algorithm corresponds to a slight variation of the multiple testing step of (Tuleai - 

III 


Malot et al. 


20141, but adapted to non necessarily symmetric distributions 


In several 


applications, neuroscientists are interested in detecting dependence periods for which 
the coincidence count is only significantly too large. In this case, one can use the re- 


'Note in particular that for a fixed W, one cannot have both < 0.5 and < 0.5 and therefore, 


if a VF is detected, it can only be because of one of the two situations, p^ < p^^^ or p^ < p'-^^ which 
cannot happen simultaneously. 
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stricted set of the p^’s. Then if the eonsidered windows are disjoint and if the spike 
trains are Poisson processes that are non necessarily stationary, the False Discovery 
Rate (FDR) of the above multiple testing procedure is mathematically provecj^to be 
controlled by q for any B >2. 

The code has been parallelized in C++ and interfaced with R. The full corresponding 
R-package is still a work in progress but actual codes are available at 

https://github.com/ybouret/neuro-stat. 


4.2 Comparison on simulations 

Two sets of simulations have been performed, the corresponding results are described 
in Table and one run of simulation of the Permutation UE method is presented in 
Figure |7J Four methods have been compared: 


the MTGAUE method of Tuleau-Malot et al. (2014) which assumes both pro¬ 


cesses to be homogeneous Poisson processes. 


the Trial-Shuffling, version C (TSC) which corresponds to the method of Pipa & 


Griin (2003), which has been programmed with the delayed coincidence count 
described above and which has not been corrected for multiplicity. 


the same as above but corrected by Benjamini and Hochberg procedure (TSC + 


BH), 


see (Tuleau-Malot et al. 


2014 1 or Table|l|for a precise definition 
^Thep(j^’s are independent random variables such that Px(p((/ < a) < aforallain [0,1] (Benjamini 


& Yekutieli 2001 Romano & Wolf 20051. 
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• the Permutation UE approaeh deseribed above. 


The permutation approaeh always guarantees an FDR less than the preseribed level of 
0.05 whereas MTGAUE does not when the homogeneous Poisson assumption fails (Ex¬ 
periment 1). The classieal trial-shuffling method (where dependenee detection occurs 
each time the p-value is less than 0.05) seems to have comparable results in terms of 
both FDR and False Non Discovery Rate (FNDR) on Experiment 1 but fails to con¬ 
trol the FDR on the most basic situation, namely purely independent processes (Ex¬ 
periment 2). Adding a Benjamini-Hochberg step of selection of p-values to the trial¬ 
shuffling makes it more robust but at the price of a much larger FNDR with respect to 
the Permutation UE method, fact which is consistent with the conservativeness shown 
in Figure 


4.3 Comparison on real data 

Behavioral procedure The data used in this theoretical article to test the dependence 
detection ability of the four methods were already partially published in previous ex¬ 


perimental studies (Riehle et al. 2000 Grammont & Riehle 2003; Riehle et al. 2006) 


and also used in (Tuleau-Malot et al. 20141. These data were collected on a 5-year-old 


male Rhesus monkey who was trained to perform a delayed multidirectional pointing 
task. The animal sat in a primate chair in front of a vertical panel on which seven 
touch-sensitive light-emitting diodes were mounted, one in the center and six placed 
equidistantly (60 degrees apart) on a circle around it. The monkey had to initiate a trial 
by touching and then holding with the left hand the central target. After a fix delay of 


500ms, the preparatory signal (PS) was presented by illuminating one of the six periph- 
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Independent 


Dependent 


spont=60 

hL^i=-600.1 


[ 0 , 0 . 001 ] 



^^‘^[0,0.005] h. . = 0 

l->J 


h. f6.1_0.005] 

|->J [0,0.005] ■ 5 


0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 


Poisson 


Hawkes 



A: Description of Experiment 1 


B: Result of the permutation method 


Figure 7: Multiple tests. BA: deseription of Experiment 1. In the Poisson part, the 
intensity of both Poisson proeesses is plotted. The injeetion eomponent eorresponds to 
the part of a shared Poisson proeess whieh is injeeted in both proeesses eorresponding 
to and X^, as explained in Figure In the Hawkes part (see ( Tuleau-Malot et al 


2014) for a eomplete deseription), formulas for the spontaneous parameters and both 
self interaetion hi^i and eross interaetion hi^j funetions are given. BB: results of 
the Permutation UE method {B = 10000, q = 0.05) performed on 191 overlapping 
windows of the form [a, a + 0.1] for a in {0, 0.01,..., 1.9} on one run of simulation for 
50 trials of Experiment 1. A blaek (resp. gray) eross is represented at the eenter of the 
window when it is deteeted by a (resp. p(y). Eaeh horizontal line eorresponds to a 
different 5 in {0.001, 0.002, ...0.04}. The blaek vertieal lines delimit the regions where 
the independenee hypothesis is not satisfied: plain for positive dependenee (i.e. where 
Qobs should be too large), and dashed for negative dependenee (i.e. where 0°^"^ should 
be too small). The dotted vertieal line separates the region of high (on the left) and low 


(on the right) dependenee in the Hawkes positive dependenee part. 
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Independ. 

Depend. 

Total 

Rejected 

V 


R 

Accepted 

U 

T 

m — R 

Total 

mo 

m — mo 

m 



Exper 

FDR 

iment 1 

FNDR 

Exper 

FDR 

iment 2 

FNDR 

MTGAUE 

0.10 

0.17 

0.04 

0 

TSC 

0.01 

0.26 

0.25 

0 

TSC -1- BH 

0 

0.32 

0 

0 

P 

0.01 

0.23 

0.02 

0 


Table 1: False Diseovery and Non Diseovery Rates. On the left hand-side, the elassieal 
table for multiple testing adapted to our dependenee framework, with a total number of 
tests m = 2K. On the right hand-side, estimated FDR and FNDR over 1000 runs, FDR 
being defined by E [(l//i?)lj:j>o] and FNDR being defined by E [(T/(m — R))lm-R>o]- 
Experiment 1 is deseribed in Figure |7} Experiment 2 eonsists in two independent ho¬ 
mogeneous Poisson proeesses of firing rate 60 Hz on [0,2]. The set of windows is as 
in Figure]^ There are 50 trials and 5 = 0.01s. MTGAUE is the method deseribed in 


(Tuleau-Malot et al. 20141 with q = 0.05. (TSC) is the trial-shuffling method with 


Monte-Carlo approximation (B = 10000) and the seleeted windows are the ones whose 
p-value are less than 0.05. (TSC-i-BH) is the same method, except that the multiplicity of 
the tests is corrected by a Benjamini-Hochberg procedure (q = 0.05). (P) corresponds 
to the Permutation UE method (B = 10000, q = 0.05). 


eral targets in green. After a delay of either 600ms or 1200ms, selected at random with 
various probability, it turned red, serving as the response signal and pointing target. 
During the first part of the delay, the probability Presp for the response signal to occur at 
(500-|-600)ms = 1.1s was 0.3. Once this moment passed without signal occurrence, the 
conditional probability for the signal to occur at (500 -f 600 -f 600)ms = 1.7s changed 
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to 1. The monkey was rewarded by a drop of juice after each correct trial. Reaction 
time (RT) was defined as the release of the central target. Movement time (MT) was 
defined as the touching of the correct peripheral target. 


Recording technique Signals recorded from up to seven microelectrodes (quartz in¬ 
sulated platinum-tungsten electrodes, impedance: 2 — 5Mf2 at lOOOHz) were amplified 
and band-pass filtered from 300Hz to lOkHz. Using a window discriminator, spikes 
from only one single neuron per electrode were then isolated. Neuronal data along with 
behavioral events (occurrences of signals and performance of the animal) were stored 
on a PC for off-line analysis with a time resolution of lOkHz. 

In the following study, only trials where the response signal (RS) occurs at 1.7s are 
considered. The expected signal (ES) corresponds to an eventually expected but not 
confirmed signal, i.e. at 1.2s. Only the pair 13 of the previous data set is considered 


here, as it was one of the main two examples already treated in (Tuleau-Malot et al 


2014). 


The results are presented in Figure]^ The (TSC-i-BH) method does not detect any¬ 
thing and is therefore not presented. The Permutation UE method detects less windows 
than both (MTGAUE) and (TSC) methods, but the detected windows are still in ade¬ 
quation with the experimental or behavioral events. The above simulation study let us 
think that the extra detections of both (MTGAUE) and (TSC) may be false positives, 
since both methods do not control the FDR as well as the Permutation UE method. 
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Figure 8: Raster plots of the pair of neurons 13. In blaek the Unitary Events where the 
eoineidenee eount is signifieantly too large for the three methods (MTGAUE, TSC and 
P) presented in Tablewith S = 0.02s and B = 10000. No interval was deteeted for 
a signifieantly too small eoineidenee eount. Signs on bottom eorresponds to behavioral 
events. The first blaek vertieal bar eorresponds to the preparatory signal (PS), the gray 
vertieal bar to the expeeted signal (ES), the seeond blaek vertieal bar to the response 
signal (RS). 
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5 Discussion 


After describing a fast algorithm to compute the delayed coincidence count, showing 
that this notion can be used in practice for any surrogate data method in place of the 
binned coincidence count, we have focused on distribution free methods to test the 
independence between two simultaneously recorded spike trains. Though they are here 
presented with the delayed coincidence count, all these distribution free methods could 
be applied to any coincidence count if desired. 

Once the coincidence count C chosen, we have first introduced an empirical quan¬ 
tity or statistic U whose distribution is centered under the independence hypothesis 


(Hq). In the spirit of (Tuleau-Malot et al. 2014| but in a distribution free manner, a 
first naive method consists in performing a Gaussian approximation of the distribution 
of U under (Hq). This method suffers from a not very sharp approximation when the 
number of trials n is small (see Figure [^. Moreover the approximation is clearly not 
valid when the observed sample does not satisfy (Hq) (see the last line of Figure]^. 

We then turned to bootstrap methods. One of the most well-known bootstrap method 
in the neuroscience literature is the trial-shuffling (Pipa & Griin 2003 [ Pipa et al. 20031. 
It is usually based on a resampling approach directly applied on the coincidence count 
itself, namely C. The other two investigated methods (full bootstrap and permutation) 
are well-known in the statistics literature for independence testing between real valued 


random vectors (Hoeffding 1952; Romano 19891 and have been recently applied to 


point processes that are modeling here simultaneously recorded spike trains (Albert 


et al. 20141. These last two methods are usually applied on centered statistics. 


One of the main message of the present work is that applying bootstrap methods to 
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non centered statisties, as C, does not lead to a eorreet approximation of the distribution 
of C under (Hq) (see the seeond line of Figure]^. This phenomenon, eombined with 
observations in the more classieal framework of the naive test on Figure leads us 
to think about a eentering defect. Once the methods are applied on correct centered 
statistics (U or U depending on the chosen bootstrap method), they all three outperform 
the naive approaeh. The approximation is better under (Hq) for small value of n (first 
line of Figure|^and first line of Figure]^ and is still aeeurate when the observed sample 
does not satisfy (Ho) (last line of Figureand second line of Figure [^. 

From an algorithmie point of view, all the corresponding bootstrap p-values are eval¬ 
uated thanks to a Monte-Carlo algorithm and a program whieh interfaees R and C++, 


thus making the running time fast and the use easy. Pipa and Grfin (Pipa et al. 2003) 
have given an exaet algorithm when the trial-shuffling is applied on the eoineidenee 
eount C direetly. It is a very elegant algorithm using the faet that C is an integer that 
can take a small number of values. Unfortunately the same gain is not really possible for 
U which is not an integer and whieh can take mueh more values. Moreover this exact 
algorithm is quite long with respeet to the Monte-Carlo algorithm when the number of 
simulations is 10000 (as used in the present work) and one ean see on the bottom left of 
Figure that the differenee between both results ( Monte-Carlo and exaet algorithms) 
is not deteetable at first glance. 

A more preeise study of the Monte-Carlo approximated bootstrap p-values shows 
that for a small number of trials, the trial-shuffling and the full bootstrap methods, even 
applied to a eorreetly eentered statistie, do not provide tests of preseribed level. On 
the eontrary, the permutation method, thanks to an adequate version of its p-values 
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(Romano & Wolf , 20051, allows for a precise control of the level. The classical trial¬ 
shuffling method based on the non centered quantity C and the naive approach also 
both lead to a precise control of the level but in a more conservative way (see the right- 
hand side of Figure]^. This is also showed by the behavior of the p-values under the 
alternative, p-values that are smaller for the permutation approach than for the other two 
methods (see the right-hand side of Figure]^. 

Finally, we decided to combine the delayed coincidence count, which is much more 


precise than the binned coincidence count (Tuleau-Malot et al. 2014; Griin et al. 1999) 


with the permutation approach, and to apply the obtained independence testing proce¬ 
dure to several windows of detection simultaneously. The final proposed method con¬ 
sists in combining the individual tests with the approach of (Benjamini & Hochberg[ 


1995) to correct for the multiplicity of the tests. Parallel programming is used to treat 
each window in an independent manner. This new algorithm named Permutation UE is 
completely distribution free. It better controls the False Discovery Rate than MTGAUE 


(Tuleau-Malot et al. 2014) or the classical trial-shuffling method applied on C (see 


Table [T] methods (MTGAUE) and (TSC)). Moreover, it does not suffer from conserva¬ 
tiveness as the trial-shuffling method applied on C, once the multiplicity of the tests is 
taken into account (see Table [TJ method (TSC-i-BH)). On real data, the results are simi¬ 
lar to existing methods (MTGAUE, TSC) except for some detections that disappear but 
that are likely to be false positive thanks to the present simulation study. 

To conclude, we introduce in this article the Permutation UE method, which is a 
Unitary Events method based on delayed coincidence count and on an evaluation of 
p-values via a distribution free Monte-Carlo approximated permutation approach. This 
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method does not suffer from any loss in synehrony deteetion as the binned eoineidenee 


eount (Griin 19961, is distribution free and in this sense upgrades (Tuleau-Malot et al. 


2014). Moreover the algorithm is fast and parallelized, and despite using a Monte-Carlo 


seheme, it ean guarantee the single tests to be of the prescribed level and the multiple 
test to control the FDR in a non asymptotic manner, therefore outperforming the trial¬ 


shuffling method (Pipa & Griin 2003 Pipa et al. 20031 in terms of both mathematical 
caution and computing time, when compared with the exact algorithm described in 
(Pipa et al. 2003[ ). Finally it is still sufficiently sensitive to detect reasonable features 
on real data sets. The only drawback is that it can only work for pairs of neurons. The 
definition of delayed coincidence count for more than two neurons has been recently 


introduced in Chevallier & Laloe (2014), but the combination of this notion with a 


bootstrap approach is still an open question. 
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